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ABSTRACT 

The traditional Schmidt density estimator has been proven to be unbiased and effective 
in a magnitude-limited sample. Previously, efforts have been made to generalize it for 
populations with non-uniform density and proper motion-limited cases. This work 
shows that the then good assumptions for a proper motion-limited sample are no 
longer sufficient to cope with modern data. Populations with larger differences in the 
kinematics as compared to the Local Standard of Rest are most severely affected. 
We show that this systematic bias can be removed by treating the discovery fraction 
inseparable from the generalized maximum volume integrand. The treatment can be 
applied to any proper motion-limited sample with good knowledge of the kinematics. 
This work demonstrates the method through application to a mock catalogue of a 
white dwarf-only solar neighbourhood for various scenarios and compared against the 
traditional treatment using a survey with Pan-STARRS-like characteristics. 

Key words: methods: statistical - proper motions - stars: kinematics and dynamics 
- stars: luminosity function, mass function - white dwarfs - solar neighbourhood. 


1 INTRODUCTION 


The use of a maximum volume as a density estimator began 
when Schmidt (19681 introduced the V/Vmaa: technique for 
analysing the luminosity function of quasars, (f)a, where V 
and Ymax are the volume enclosed by the object and the 
volume in which the object can be found at the given sur¬ 
vey limits respectively and <j)s is the number density. The 
estimator can be used as a means of testing the complete¬ 


ness of a sample simultaneously. Felten (19761 showed that 
(j)s is unbiased and superior to the older “classical” estima¬ 
tor (N/V) when the magnitude or luminosity bins are not 
small. He further described the procedure of combining the 
density estimator from two non-overlapping areas. 

Information from a single survey is limited due to the 
relatively small number of objects. When several complete 
samples are combined, more fundamental parameters of the 
population of objects can be determined, and with smaller 
associated uncertainties. In the light of this problem, |Avni| 
|fc Bahcall| ( |1980[ ) investigated different ways of combin¬ 
ing catalogues, namely the incoherent region-independent 
method, the incoherent domain-independent method and 
the coherent method. All of these are superior over the orig¬ 
inal Schmidt method, with the coherent method being most 
accurate. When the effects of a space-density gradient were 
corrected for ([Stobie, Ishida, fc Peacock|1989[ [Tinney, Reid,| 
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|fc Mould||1993| |Lam fc Hambly||2015| ) this method was ex¬ 
tended to estimate stellar density where the density profile of 
the Galaxy varies significantly along different lines of sight. 
Because of the small distances probed, only the scale height 
effects are considered while the scale length is assumed to 
be constant. 


In order to consider a sample of proper motion objects. 


Schmidt (1975 I extended his estimator to cope with both 


photometric and proper motion detection limits. The new 
estimator considers the tangential velocity as an intrinsic 
property of an object such that it can be kept as a constant. 
Then, the distance limits can be found easily by applying 
the upper and lower proper motion limits of the survey to a 
simple relation between tangential velocity, proper motion 
and distance: Vtan = 4.74/rT> kms“^, where fj. is the proper 
motion in arcseconds per year and D is the distance to the 
object in parsecs. 

Cool white dwarfs (WDs) and subdwarfs (sds) have sim¬ 
ilar optical colours to the main sequence stars (MSS) while 
those of brown dwarfs (BDs) are similar to the giants. There¬ 
fore, it is difficult to distinguish them in colour-colour space. 
Due to their small radii, WDs, sds and BDs are located far 
from MSSs and giants in the HR diagram. However, when 
objects are only detected in a few broadband filters, it is im¬ 
possible to classify them reliably which would lead to poor 
object selections and distance estimates. To overcome this 
problem, it is common to use reduced proper motion (RPM) 
as a crude estimate of absolute magnitude to separate sam- 
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pies of subluminous objects from higher luminosity contam¬ 
inants. In order to obtain a clean sample of WDs (for ex¬ 
ample an extreme subdwarf would easily be confused with 
a WDs with low tangential velocity) a lower tangential ve¬ 
locity has to be applied to remove the ambiguous objects. 
This procedure introduces an incompleteness which has to 
be corrected for. This problem was identified by |Bahcall fc| 
Casertano (19861 and Evans (19921 separately. The former 


Table 1. Physical properties of the Galaxy used in the Monte 
Carlo simulation. 


adopted a Monte Carlo (MC) simulation approach to cor¬ 
rect for the incompleteness, while the latter was done ana¬ 
lytically. However, the two methods evolved separately. In 


the simulation front, Liebert et al. (1999) constructed some 


simulations based on different Galactic models to study the 
incompleteness due to proper motion selection after some 
strong arguments ( Oswalt fc Smith] 1995 Oswalt et al.[199'6 1 
pointing towards an incomplete LHS catalogue used in ear¬ 
lier studies (eg. Winget et al. 1987 Liebert, Dahn, & Monet 


Parameter 

Thin disc 

Thick disc 

Stellar Halo 

(U}/km 

-8.62“ 

-ll.O'' 

-26.0'' 

(V}/km 

-20.04“ 

-42.0^* 

-199.0'' 

(W)/km s-i 

-7.10“ 

-12.0'* 

-12.0'' 

ajj/km s~^ 

32.4“ 

50.0'* 

141.0'' 

(Ty/km 

23.0“ 

56.0'* 

106.0'' 

iTvv/km 

18.1“ 

34.0'' 

94.0'' 

H/pc 

250'’ 

780“ 

00 

n/pc“® 

0.00310“ 

0.00064“ 

0.00019“ 


a Fuchs, Jahreifi, &: Flynn|2009| 


1988 ). This correction, known as the discovery fraction, x, 
was then applied by Harris et al. ( 2006] her eafter H06). On 
the other hand, Knox, Hawkins, & Hambly (19991 used the 


b Mendez’^r'^uzmanlT^SI 

c RHll 

d IChiba &: Beersl2000l 

e IGirard et al. 120061 

the halo has a star burst of duration 1 Gyr at r = 12.5 Gyr. 


analytical approach to correct for the incompleteness. In¬ 
stead of calculating the discovery fractions from integrating 


The initial mass function has an exponent of -2.3 (Kroupa 


20011, and the initial-final mass function(IFMF) follows the 


over the density profile, Digby et al. (20031 arrived at the dis¬ 


Kalirai et al. (20091 


covery fractions by integrating over the Schwarzschild distri¬ 


bution functions. Rowell & Hambly (2011 hereafter RHll) 


further generalized the technique to cope with an all sky 
survey as opposed to the individual fields of view employed 
in earlier works. 

This work studies the discovery fraction in detail and 
shows that the discovery fraction has to be incorporated into 
the volume integral in order to arrive at a correct density 
estimation. In the next Section, we discuss the method of 
the simulation of the solar neighbourhood. In Section 3 we 


I O.lOlrui -|-0.463, 0.5Mq < ^ d.OM© 

^ |o.047mi -k 0.679, 4 .OM 0 < mi ^ 7.0Mq. 

The MS life time has to be added in order to calculate 
the cooling time, and hence the magnitude of a WD. We 
have adopted the stellar evolution tracks from the Padova 
group (PARSEC; Bressan et al. 2012[) with a metallicity 
of Z = 0.019 and Y = 0.30 (Girardi et al. 2000). To¬ 


gether with the pure hydrogen WD (DA) atmosphere cool- 


discuss different maximum volume estimators and discov- 

synthetic colour^ 

(Holberg & Bergeron 

2006 Kowalski & 

ery fractions and how their shortcomings can be removed 
by a new approach. This new method is then applied to 

Saumon||2006| 

Tremblay, Bergeron, & G 

ianninas||2011| and 

Bergeron et al. 

20111, a theoretical luminosity function (LF) 


we choose survey parameters typical of the state-of-the- 


art Pan-STARRS optical sky survey (Hambly et al. 2013 


and references therein). In the last Section, we discuss the 
possible extension of the method and conclude this work. 


2 POPULATION SYNTHESIS 

Monte Carlo simulations are used to produce snapshots of 
WD-only solar neighbourhoods which carry six dimensional 
phase space information. The volume probed in this work 
is assumed to be small such that the simulation is done in 
a Cartesian space, instead of a plane polar system centred 
at the Galactic Centre (GC). The Galaxy is further assumed 
to have three distinct kinematic components: a thin disc, a 
thick disc and a stellar halo, all of which have no density vari¬ 
ations along the co-planar direction of the Galactic plane. 
All vertical structures follow exponential profiles, with scale 
height H. The velocity components, U, V and W, of each ob¬ 
ject are drawn from the Gaussian distributions constructed 
from the measured means and standard deviations of the 
three sets of kinematics that describes the three populations 
in the solar neighbourhood. The thin and thick disc popula¬ 
tions are assigned with constant star formation rates since 
look back time, r = 8 Gyr and r = 10 Gyr respectively while 


is produced which can then be used as the PDF in the 
Monte Carlo simulation. The normalisations of the PDFs 
are adopted from the WD densities found in RHll. The 
volume in which objects are distributed is limited to half a 
magnitude deeper than the maximum distance at which the 
survey can probe given its brightness. The half magnitude 
is to allow for random fluctuations near the detection limit 
after noise is added. The input parameters are assumed to 
be invariant with time and are summarized in Table 1. 

From the true distance and true bolometric magnitude 
drawn from the PDF, the true apparent magnitudes in the 
Pan-STARRS gpi, rpi, ipi, zpi and ypi filters are assigned 
( Tonry et al.|2012 Schlafly et al. |20T2 Magnier et al.|2(5l3 |. 
The uncertainties in those filters, am.i , are assumed to scale 
exponentially with magnitude and are described by 


= ai X exp 


- 15 . 0 ) 


-\-bi 


( 2 ) 


where ai and bi are constants measured from the first Pan- 
STARRS survey (PSl; magnitude subscript ‘PI’) at Pro¬ 
cessing Version (PV) 1.1 and rui is the magnitude in filter 
i (Fig. 1). Realistic dispersion is added to the uncertainties 
by resampling with a Gaussian distribution with stan¬ 
dard deviations of 0.1 x am centred at the noiseless am ■ 


http://www.astro.umontreal.ca/~bergeron/CoolingModels/ 
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Table 2. Parameters for the noise model. 


Filter or Proper Motion 

ail mag 

hi! mag 

g 

0.000125 

0.00065 

r 

0.000120 

0.00075 

i 

0.000175 

0.00065 

z 

0.000325 

0.00089 

y 

0.000750 

0.00120 

p 

0.000300 

0.00050 



12 14 16 18 20 22 

MagnlTide / mag 


0.10 


0.08 


0.06 


0.04 


t 0.02 


0.00 


0.10 


J 0.00 


0.06 


0 114 


0.02 


Figure 1. Top: the magnitudes and the associated uncertainties 
of point sources with 5cr proper motions in the PVl.l of Pan- 
STARRS, from right to left gpi, rpj, ipi, zpj and ypi filters. 
Each successive filter is offset by one magnitude. Bottom: the 
magnitudes and uncertainties distribution of WDs reproduced in 
the simulation. 


The magnitudes in each filter are then drawn from a Gaus¬ 
sian distribution with a standard deviation Um; • The proper 
motion uncertainty is based on the rpi magnitude. 


3 GENERALIZED MAXIMUM VOLUME 
DENSITY ESTIMATOR 

The classical estimator il? = A/V for a volume-limited sam¬ 
ple is of little practical use for analysing small numbers of 


objects or strongly localized and kinematically biased groups 
of stars that are selected by proper motion and apparent 
magnitude. For these samples, the 1/Vmaa: method is gen¬ 


erally regarded as a superior estimator of the LF (Felten 


19761. The contribution of each object to the LF is weighted 


by the inverse of the maximum volume in which an object 
could be observed by the survey. For example, for a given 
bin of objects with index k, the space density is the sum of 
all 1/V max 


-h) 


Nk 

-E 

i=l 


v„ 


(3) 


for Nk objects in the k‘" bin. The uncertainty of each star’s 
contribution is assumed to follow Poisson statistics. The sum 
of all errors in quadrature within a luminosity bin is there¬ 
fore. 


^ k 


'Nk 

E 

i = l 




(4) 


The traditional l/Vmai technique assumes that objects are 
uniformly distributed in space. However, in reality, stars in 
the solar neighbourhood are concentrated in the plane of 
the disc. The effects of space-density gradient can be cor¬ 
rected by assuming a density law and defining a maximum 
generalized volume Vgen (|Stobie, Ishida, fc Peacock||1989| 


Tinney, Reid, & Mould 19931 which is calculated by inte- 


and 


(5) 


grating the appropriate stellar density profile p/p© along 
the line of sight between the minimum distance, dmi 
maximum distance, dmax- This leads to the integral 

Vge,..S89 = n / ^ r" dr 

where Q. is the size of the solid angle of the survey. To min¬ 
imise the contamination from extreme subdwarfs scattered 
into the WD regime in RPM-colour space, a lower tangen¬ 
tial velocity limit, Vtan,lower, is applied to remove most of 
the contaminants. Traditionally, the discovery fraction, Xu 
which is the fraction of objects with tangential velocities 
larger than the lower tangential velocity limit, is only Galac¬ 
tic model and survey-footprint dependent ([Bahcall fc Caser-| 


tano||19^ Liebert et al. 1999 Harris et al.|2006 l such that 

the maximum volume density estimator can be written in 
the form 


V. 


[ 

gen,H 06 ~ \v(Vtan,lower') H / 

Jd„ 


Pjr. 

Po 


LlLf fij. 


where the distance limits are derived from both photometric 
and proper motion limits of the survey by calculating 


dr, 


= d X max 


10^ 


P 

fd'max 


(7) 


dmax = dx min 


10 


5 


fd^min 


( 8 ) 


where mmin,i, rnmax,i and Mi are the photometric limits 
and the absolute magnitudes of the object in filter i respec¬ 
tively. The proper motion terms, iJ,/^max and are 

rationalized by assuming an object would carry the same 
tangential velocity if it were placed closer to or farther from 
the observer (analogous to the absolute magnitude) and/or 
in an arbitrary line of sight. 
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3.1 Attempt to Modify the Discovery Fraction 


RHll extended the Xv to include a directional dependence in 
order to account for the varying survey properties and stel¬ 
lar tangential velocity distribution across the sky. In RHll, 
the Schwarzschild distribution function is used instead to 
calculate the tangential velocity distribution, P{vtan), ana¬ 
lytically. The discovery fraction can be found by projecting 
the velocity ellipsoid onto the tangent plane of observation 
( Murray! 19831. Thus it allows one to arrive at a precise 
x(a, d) without taking the average properties over a large 
area. Therefore, the volume integral can be modified to 


V, 


gen,RHll 




'^tan,lower, 


f 


P^r^dr (9) 
Pq 


where i denotes each sky cell covered by a Schmidt survey 
field employed in the production of the catalogue used by 
RHll (hereafter the generalized method). This should, in 
theory, have taken into account all the small scale variations 
which a positional-independent would not be able to deal 
with. 


3.2 Problems in this framework 

We have identified a bias present in the whole family of max¬ 
imum volume methods when calculating distance limits by 
holding the tangential velocity constant along a single line of 
sight (Schmidt 19751; and the consequence to the the discov¬ 
ery fraction applied in HOG and RHll when the tangential 
velocity is further held constant across the sky. These can 
be described as follows: 


Constant tangential velocity along line of sight 

The kinematics of an object is a property of the Galaxy. An 
object at a given magnitude at any given distance from the 
observer should not carry the same tangential velocity at a 
different line of sight distance when tested for observabil¬ 
ity. This assumption is only good over small field of views 
and small range of line of sight distances, while an accept¬ 
able size of the smallness is very difficult to access if not 
unquantifiable. 


Constant tangential velocity across the sky 

Consider a spatially uniform population like the stellar halo, 
and an all-sky survey where the proper motion limits are the 
same over the whole sky. The tangential velocity distribution 
varies along different lines of sight due to the solar motion. 
For stars at a given magnitude, a different fraction of the 
population will pass the proper motion limits along differ¬ 
ent lines of sight due to the differences in the tangential ve¬ 
locity distributions. In the most extreme cases, on average, 
a halo WD observed in the direction of the Anti-Galactic 
Center (AGC) would appear to have a large velocity due to 
reflex motion imparted by the Sun in its orbit within the 
Galaxy. However, if one is observed in the direction to the 
solar apex instead, the motion would be much smaller on 
average. The consequence is that when proper motion limits 
and tangential velocity limits are applied, different numbers 
of stars would be detected in different regions of sky even for 



Figure 2. This plot illustrates how tangential velocity and 
proper motion limits behave in the y — D space. The solid lines 
are the contours of the tangential velocities from 20 to 200kms“^ 
in steps of 20kms“^. The photometric distance limits mark the 
range of distances in which an object can be placed and stay 
within the detection limits. The dot-dashed lines are the distance 
limits calculated from Equations and which is by fixing the 
tangential velocity of an object such that a proper motion limit 
correspond to a fixed distance limit. The dotted lines are the 
distance limits calculated from Equations |13| and |14| where the 
distance limits are not functions of the kinematics. The dashed 
lines are the upper and lower proper motion limits; in the lat¬ 
ter case, this limit rises as astrometric errors rise for increasingly 
faint flux levels as distance increases. 


identical survey limits and spatial density. That is the moti¬ 
vation for RHll to consider different tangential velocity dis¬ 
tributions along different lines of sight. The method would 
be completely correct in the framework where the tangential 
velocity is an intrinsic property of an object. However, in the 
case where tangential velocity is not constant, there would 
be a mismatch between the parameter space which the dis¬ 
covery fraction and the maximum volume explore (see Sec¬ 
tion where we will demonstrate how the new approach 
can solve both problems). 


3.3 A Closer Look at the Discovery Fraction 

Proper motions, tangential velocities and distances are re¬ 
lated by a simple equation 

uton ~ 4.74/rD. (10) 

All the survey limits can be shown in the proper motion- 
distance space (Fig. [^. A valid approach should not im¬ 
pose any assumptions that restrict an object in this two 
dimensional space. The discovery fraction is equal to the 
weighted area restricted by the survey limits, where the 
weight map (Fig. is generated from dividing the tangen¬ 
tial velocity distribution by the distances, 

W{fi,D) = (11) 

4.74 X D 

In the generalized method, Vtan is fixed and is thus separable 
from fi, which follow the lines of constant tangential veloc¬ 
ity at Vtan = 20 - 200 kms“^ in 20kms“^ interval in Fig.[^ 
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Figure 3. The area in light grey (which also includes the hatched 
and dark grey areas) is the parameter space where the discovery 
fraction is calculated when the proper motion limits and distance 
limits are not considered (e.g. H06 and RHll). The dark grey 
area corresponds to the generalized volume where the distance 
limits (Eq. and are considered, it is clear that there are 
inconsistencies between the two parameter space. The hatched 
area is the parameter space where both the discovery fraction and 
the volume are integrated over in the modified method (Please 
note that in this case the distance limits are found from Eq. |13| 
and |14| l. The weighted area gives the discovery fraction in the new 
method (see Figure|^for the weight maps). A treatment without 
considering the effects of the proper motion limits and distance 
limits on the parameter space in which the discovery fraction 
is integrated over would always overestimates the completeness, 
which translates to an underestimation in the number density. 



Distance / pc 


However, the two variables are related through distance as 
shown in Equation |10[ This implies a selection criterion in 
one parameter would lead to a selection effect in the other. 
The generalized method only deals with the tangential veloc¬ 
ity limits but ignores the consequential effect of the change 
in the proper motion limits bounding the discovery fractions. 
A similar effect also appears in the treatment of the distance 
limits. The volume integral is bounded by the maximum and 
minimum distances in which an object can be found but the 
discovery fraction in RHll includes everything above the 
line vtan = 20kms“^ as represented by the light grey area 
in Fig.|^ The true discovery fraction should be bounded by 
the same limits as applied to the volume integral (hatched 
area). Thus, Xv is always overestimated which translates to 
an underestimation in the luminosity function. This effect is 
stronger for: 

i) a survey with small upper proper motion limit because 
the lower this limit, the larger the overestimation of the dis¬ 
covery fraction. This can be seen from Fig.j^where a smaller 
upper proper motion limit would lead to a smaller cross- 
hatched area, which means the discovery fraction would be 
overestimated. 

ii) a population with large differences in the kinematics 
compared to the observer because objects tend to have large 
proper motions. This shifts the region with the highest prob¬ 
ability in W{fi,D) (Fig.|^ to larger proper motions. In the 


Figure 4. The top panel shows the probability distribution as a 
function of proper motion and distance of a volume limited thin 
disc population in the direction of (a, S) = (0.0, 0.0) over-plotted 
with the proper motions, tangential velocities and distance limits. 
The discovery fraction is the weighted area within the specified 
limits where the probability distribution function is used as the 
weights. The lower panel shows those of a halo population. The 
heaviest part of the thin disc weight function is excluded by the 
proper motion limit, but much more of the heavily weighted re¬ 
gion of the stellar halo is within the allowed parameter space, so 
the tangential velocity limit at ~ 20kms“^ has a much smaller 
effect in the thin disc than in the halo. 


case of the Galaxy, the stellar halo is the most susceptible 
to this effect. 

iii) intrinsically faint objects that carry small maximum ob¬ 
servable distance. It is clear from Fig.[^ that the relative 
W{j-L, D) at low proper motion increases with distance when 
the P(vtan) / D peaks at smaller proper motion as distance 
increases. This affects the discovery fraction more severely 
when the lower proper motion limit is large since more 
heavy-weighted area would be included in, for example, both 
the H06 and RHll methods. 


As an example of a current survey which can be employed 
in WDLF studies (Hambly et al. 20131, the Pan-STARRS 
upper proper motion limit is 0(1)” yr“^ which means the 
first problem may affect the analysis of WDLF without 
proper treatment. Halo WDs have higher velocities and are 
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older (i.e. fainter) than those of the discs and hence the effect 
on the halo population would be much larger than that in 
the discs. As it is of great interest to probe the halo WDLF 
at such low luminosity to explore the possible scenarios of 
the star formation history of the Galaxy, it is necessary to 
correct this bias. 


3.4 A New Approach 


In order to compute the discovery fraction properly, the pa¬ 
rameter space in which an object could be observed by the 
survey has to be identical to that used in the discovery frac¬ 
tion integral. For each step of the numerical integration, it 
is necessary to calculate the instantaneous discovery frac¬ 
tion which is limited by the upper and lower proper motion 
limits, as well as the tangential velocity limits. It is worth 
mentioning that in the effective volume method in RHll a 
similar approach was adopted that would have corrected for 
the bias noted in the last Section although the bias was not 
explicitly identified and discussed in that work. Instead of 
dealing in an object by object basis, their correction was ap¬ 
plied statistically. The strength of that method is that the 
WDLFs of the three components could be untangled. How¬ 
ever, binning objects by their bolometric magnitudes before 
membership association would lead to a loss of information. 
Furthermore, it loses the generality so that it cannot be ap¬ 
plied to other luminosity estimators. In order to keep the 
analysis in an object by object basis, one should consider 
the modified volume integral. It is computed by integrating 
the physical survey volume along the line of sight, and con¬ 
sidering at each distance step both the stellar density profile 
and the fraction of objects that pass the tangential velocity 
limits implied by the proper motion limits. The total gen¬ 
eralized survey volume between dmin and dmax is therefore 
calculated by 





rb(r) 

1 P{VT)dVT 

■'dmin ^0 

<i(r) 


dr, ( 12 ) 


where the distance limits are solely determined by the pho¬ 
tometric limits of the survey 


10 " 


i-M,) 


(13) 


dmax = d X min 



(14) 


P{vt) is the tangential velocity distribution, p{r) corrects 
for a non-uniform population density profile (i.e. disc pop¬ 
ulations), and Q is the survey footprint area in steradians. 
Note that P{vt) and p(r) depend on the line of sight, so 
this model only holds for small fields. The limits on the tan¬ 
gential velocity integral are 

a(r) = max(uinin,4.74/iininr) (15) 

&(r) = min(umax,4.74/rniaxJ'), (16) 

where Umin and Umax are any fixed tangential velocity limits 
applied to reduce contamination from other stellar popula¬ 
tions, and 4.74/rminr and 4.74/rmaxr are the tangential ve¬ 
locity limits at distance r arising from the proper motion 
limits. The appropriate limits on the integral are found by 
considering both of these effects. 


The new volume integral has the distance limits de¬ 
coupled from the kinematics, which are completely ab¬ 
sorbed into the discovery fraction. The decoupling simul¬ 
taneously means that regardless of how the kinematic be¬ 
haviour changes with respect to the direction of observation, 
the entire p-D-vtan parameter space in which an object can 
be found is explored. When the discovery fraction has a de¬ 
pendence on distance, it cannot be separated from the in¬ 
tegrand. This means that the discovery fraction varies from 
object to object in any given direction, so it is sufficiently 
general to take a more realistic form of velocity distribu¬ 
tion. This in turn allows for a distance dependent tangential 
velocity distribution that describes the Galaxy more realis¬ 
tically. 


Survey volume generalized for kinematic selection 

It is instructive to examine how the survey volume as a 
function of distance is changed when the full effects of kine¬ 
matic selection are taken into account, i.e. considering that 
proper motion limits result in implicit tangential velocity 
limits that vary as a function of distance. The generalized 
1/Vmax method does not consider this, and computes a con¬ 
stant discovery fraction only from any fixed, external tan¬ 
gential velocity limits that are applied e.g. to reduce con¬ 
tamination from certain types of object. 

The differential survey volume as a function of dis¬ 
tance, for the stellar halo and the thin disc, is presented 
in Fig. 1^ These plots are computed for a line of sight to¬ 
wards the NGP, in order to exaggerate the effect of the non- 
uniform population density profile in the disc case, which 
falls off most rapidly in directions perpendicular to the 
plane. Halo and thin disc kinematics follow those in Table 
1, with lower/upper proper motion limits of 0.1/1.0”yr“^ 
and a survey solid angle of 0.01 sr. 

In the halo case, a fixed lower tangential velocity limit of 
Umin=200kms“^ has been applied, which is the usual way 
of reducing contamination from disc WDs. Umax has been 
left unconstrained; this corresponds to a constant discovery 
fraction of ~0.67, which is the value used in the generalized 
1/Vmax method leading to a generalized volume that is a 
constant fraction of the physical volume. 

In the disc case, a fixed lower limit of Umm=30kms“^ 
has been applied, which is used to reduce contamination 
from high velocity subdwarfs. Umax is again unconstrained, 
leading to a constant discovery fraction of ~0.66 used by the 
generalized 1/Vmax method, although note that in the disc 
case the generalized volume diverges from the true physical 
volume due to the non-uniform population density profile. 

In both cases, there is a range of intermediate distances 
over which the generalized and modified 1/Vmax methods 
give the same result for the survey volume. This corresponds 
to the range over which the fixed external tangential veloc¬ 
ity limits are active. Significant differences between the two 
methods arise at large distances, where the lower tangential 
velocity limit implied by the lower proper motion limit ex¬ 
ceeds the fixed external limit. It is at these distances that 
the generalized 1/Vmax method, which fails to consider this 
effect, overestimates the generalized survey volume and thus 
underestimates the number density. A similar effect arises at 
small distances (inset), where the upper tangential velocity 
limit implied by the upper proper motion limit is greatly re- 
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Figure 6. Top: WDLFs for a population with halo kinematics 
with the measured halo WD density. The two WDLFs are derived 
from the new approach (solid) and the RHll method (dashed) 
overplotted with the input true WDLF (grey). The lines above 
and below the WDLFs are the 1-sigma upper and lower limits. 
Middle: the deviations of the WDLF as a function of magnitude. 
The RHll solution departs from the true LF at ~12.0mag while 
the new approach at ~15.0mag. Bottom: the (V/Ymax) as a 
function of magnitude from the two methods. The uncertainties in 
(V/Ymax) is i'll® (Y/Ymax)s are within the statistical 

fluctuations up to 16 mag. 


Figure 5. generalized survey voiume as a function of line of 
sight distance, for the steiiar halo (top) and thin disc (bottom). 
The lines represent the true physical volume (solid), the gener- 
aiized volume (dashed), and the modified voiume (dotted). The 
functions diverge at large distances where 4.74/rniinr > iimin (see 
text) and the discovery fraction used by the generalized 1/Vmax 
method breaks down. A similar effect operates at small dis¬ 
tances (inset), where the upper proper motion limit results in 
a small upper tangential velocity limit that excludes a large frac¬ 
tion of the population. 


duced, causing most of the population to be excluded. The 
generalized 1/Vmax method, which considers only the fixed 
threshold Umax, fails to account for this effect and again over¬ 
estimates the survey volume. 


4 APPLICATION TO WHITE DWARF 
LUMINOSITY FUNCTIONS 

We have used the simulation outlined in Section to gen¬ 
erate mock Pan-STARRS catalogues to compare the differ¬ 
ences in applying the RHll generalized volume and the new 
modified volume described in this work. The bright limits in 
all filters are set at 15.00mag, while the faint limits in gpi, 
rpi, ipi, zpi and ypi filters are at 21.50mag, 21.00mag, 
20.50 mag, 20.00 mag and 20.00 mag respectively. The lower 
proper motion limit is defined as five times the proper mo¬ 


tion uncertainty at the given magnitude as defined in Equa¬ 
tion The upper proper motion limit is set at 0.3” yr“^ 
unless specified otherwise. Photometric parallaxes are found 
by fitting the magnitudes to the WD synthetic atmosphere 
models at fixed surface gravity logg = 8.0. Tests have shown 
that there exists local minima for effective temperature be¬ 
low 6000K so a Markov-Chain Monte Carlo method is used 
for minimisation. T he implementation adopted th e Python 
programme emce^ ( [Foreman-Mackey et al.||2013 1. A lower 
tangential velocity limit is set at 20kms“^. 


4.1 Stellar halo 

As discussed in Section [3.3| faint objects and objects from 
a population with large difference in the kinematics from 
the observers are most susceptible to underestimation of 
the number density when using the generalized volume tech¬ 
nique. Thus, the largest systematic errors are expected to be 
found in the faint end of the stellar halo WDLF. In the top 
panel of Fig|^ the differences in the WDLF constructed by 
the two methods are shown. The two LFs agree with each 
other up to Mf,o; ~12.0mag. Beyond that, the generalized 
method consistently underestimates the number density and 


^ http://dan.iel.fm/emcee/ 
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Figure 7. Same as Fig|6] except the WD density is artificially 
inflated by 100 times to reduce random noise. With the increase 
number density, the WDLF departs a magnitude fainter in both 
cases. The (V/Ymax) of the old method is still fluctuating about 
0.5 even when the WDLF departs further and further away from 
the true LF. In the new method, the (V/Ymax) is consistent 
within statistical uncertainties. 


the deviation increases as the objects get fainter. The maxi¬ 
mum difference is more than 1.0 dex. The Pan-STARRS fil¬ 
ter gpi reaches ~400nm at the blue edge, so it is expected 
that the photometric parallax solutions become unreliable 
at the bright end, while with small numbers of objects in 
the faint end the statistical noise is significant. Most of the 
{V/Ymax)s in both cases are within Icr from 0.5, which is a 
necessary condition for an unbiased sample. 


4.2 Stellar Halo with 100 times the Observed 
Density 

The stellar halo has a very low number density and it is 
well known that maximum volume estimators are prone 
to systematic bias with small number statistics. Thus, we 
have generated a stellar halo with 100 times the observed 
density to reduce such uncertainties. The WDLFs pro¬ 
duced are more easily compared when the systematic un¬ 
certainties are much smaller and the corresponding plots 
are shown in Fig. The two LFs agree with each other 
up to Mboi ~ 12.0 mag again. With 100 times the density, 
faint objects down to 17.0 mag (Tefi ~ 3000 K) can be gen¬ 
erated easily. The LF produced with the new method agrees 
with the input LF to luminosities down to ~16.0 while the 
RHll method consistently underestimates the density, as 
predicted in Section [3.3[ The (V/Ymax) shows an interesting 
behaviour - the distribution is fluctuating about 0.5 regard¬ 
less of the bias in the WDLF. It implies that (V /Ymax) ~ 0.5 
is a necessary but not sufficient condition for both bias and 
completeness. 


4.3 Different upper proper motion limits 


As mentioned in Section [3.3| a small upper proper motion 
limit would lead to a strong bias. This is illustrated in Fig|^ 
as the upper proper motion limit increases, the old and new 
approaches converge. However, when the limit decreases, 
the underestimation in the number density increases when 
adopting the old method. From the bottom panel, even at 
an upper limit of 0.5” yr“^, the generalized method fails to 
recover the WDLF beyond ~15 mag. With the new method, 
the number density estimation is recovered under any re¬ 
stricted proper motion selection within the statistical noise. 
In the modern and future surveys with the rapid photomet¬ 
ric systems, the upper proper motion limit would be in the 
order of arcseconds so this effect would be small. However, 
when the pairing of the bright objects includes older pho¬ 
tographic plate data, the upper proper motion limit would 
be severel y hampered due to a large maximum epoch differ¬ 
ence (e.g. Digby et al.|2003 has an upper limit of 0.18” yr“^ 
from pairing between POSS-I, POSS-H and SDSS). Another 
case is when samples at different velocity ranges are ana¬ 
lyzed separately (e.g. the effective volume method in RHll 
where the discovery fraction was treated correctly) the re¬ 
stricted velocity range would amplify the shortcoming of the 
old method. 


4.4 Thin Disc, Thick disc and Stellar Halo 

In the study of the Galactic WDLFs it is very difficult to as¬ 
sign population membership to individual objects. The con¬ 
sensus is that thin disc objects dominate the solar neigh¬ 
bourhood so that when studying the LF for the thin disc it 
is possible to assume a thin disc characteristic for all objects. 
This does not happen to be a good assumption as shown in 
the top panel of Fig.[^where for the purposes of this exercise 
we apply only the new method. The solid line is the WDLF 
constructed from all stars regardless of the population. It is 
overestimated by ~0.2 dex at all magnitudes. An overesti¬ 
mation is expected but a consistent overdensity of 0.2 dex 
is not negligible. Due to the spatial density correction, the 
maximum volume of an object integrated over a disc profile 
has smaller volume than a halo profile. While objects are 
weighted by the inverse maximum volume, the apparently 
negligible contribution from the older populations would by 
amplified by density correction such that each contaminat¬ 
ing object would have a volume underestimated by tens of 
percentage points. The discovery fraction for a disc popula¬ 
tion is always smaller than that of the halo, so this amplifies 
the XjYraax by another few tens of percentage points. The 
YjYmax distributions are consistently at around 0.55 as a 
consequence of the halo objects contribution. The ratio of 
the volume element between an exponential disc density pro¬ 
file and a uniform halo increases with distance. Therefore, 
for a group of uniformly distributed objects being assigned 
to follow an exponential profile in the integrals, the resulting 
ly jYmax) is expected to be larger than 0.5. The dashed line 
is constructed with an addition of a fixed upper tangential 
velocity limit at 80kms“^ which would eliminate most of 
the halo objects, but not the ones from the thick disc. This 
reduces the contamination down to roughly 0.1 dex at all 
magnitudes, with a smaller effect on the {Y/Ymax) distri¬ 
bution. The dotted line is computed using only thin disc ob- 
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Figure 8. Top: WFLD constructed with different upper proper 
motion limits at 100 times the measured density with the new 
modified volume method. From top to bottom, the upper proper 
motion limit is set at 10.0, 1.0, 0.5, 0.4, 0.3, 0.2 and 0.1” yr“^. 
Bottom: same procedure repeated with the generalized volume. In 
the two sets with the smallest limits, dashed lines are used to join 
up the LFs at the faintest bin where the area under the dashed line 
should not be included in determining the total number density. 


jects from that simulation for comparison. Everything agrees 
to within statistical uncertainties at all magnitudes in the 
thin disc-only analysis. Since the discovery fraction requires 
the kinematics of the populations, a good membership asso¬ 
ciation is very important. 

The aforementioned RHll effective volume method in 
Section 13.41 was not dealt with in detail because of some 
fundamental differences in the method. That method can 
untangle the three components of the Galaxy. However, the 
binning of objects at an early stage means that it is not di¬ 
rectly compatible with the framework in this work which, 
instead, bins objects in the final step. In dealing with ob¬ 
servational data without the labels “thin disc”, “thick disc” 
and “halo” tagged on the objects, the best one can do is 
to perform the membership assignment statistically. The ef¬ 
fective volume method would provide a quick and easy way 
to do the job. However, if one would like to retain as much 
information as the maximum volume family would normally 
allow, the best way to do it is to find the probability of 



Figure 9. Top: WDLFs of the Galaxy constructed by assum¬ 
ing negligible halo contamination (solid), with an additional fixed 
upper tangential velocity limit at 80 kms“^ (dashed) and by pick¬ 
ing up only objects tagged as “thin disc object” from the simula¬ 
tion (dotted). Middle: the differences in the WDLF from the true 
LF as a function of magnitude. The overdensities are expected 
from the wrong association of kinematics model with the objects. 
Bottom: the (V/Vmax) distributions as a function of magnitude. 
{V/Vmax)s are consistently larger than 0.5 implying objects tend 
to be found at larger distances than on average. This is expected 
from the objects with a uniform spatial density being associated 
with ones that have an exponential profile. The larger discovery 
fraction of the thin disc with respect to the halo amplifies the the 
effect further (see text for detailed explanation). 


each object belonging to each of the components based on 
their observed and intrinsic properties, e.g. tangential veloc¬ 
ities metallicities, and/or distance from the Galactic Plane. 
Population membership assignment is beyond the scope of 
this paper but if it were to be done in a maximum likeli¬ 
hood approach, readers are reminded that distributions are 
usually approximated as Gaussian functions. In the case of 
tangential velocity, the PDF in fact follows a Schwarzschild 
distribution in which case a Gaussian distribution cannot 
approximate the tail at high velocity well. It is important to 
assess whether this would lead to a bias in assigning objects 
to higher velocity populations. However, if one were to split 
the tangential velocity into {l,b) components, the distribu¬ 
tions of the component velocities should be better approxi¬ 
mated by Gaussian functions. As for the physical position, 
both the current distance from the plane and the maximum 
distance in which an object can reach based on the current 
kinematics, Zmax, would be good parameters to test. Both 
suggested parameters have used the important piece of in¬ 
formation that a halo object is much more likely to carry a 
larger velocity perpendicular to the plane than those of the 
discs. 
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Figure 10. Top: WDLFs for the stellar halo with different input 
Galactic models. The solid line is the original halo WDLF iden¬ 
tical to the one in Fig[7| the dashed line is the one that follows a 
simplified Besangon model but then analyzed with the original set 
of kinematics. Middle and bottom panels are the same as before. 


4.5 Different Galactic Models 


The accuracy of all LF estimators is sensitive to the as¬ 
sumed population kinematics and density profiles. There¬ 
fore it is important to investigate the effects of apply¬ 
ing a different kinematic model in the analysis of a pop¬ 
ulation simulated with another model. To demonstrate 
the effect, we have adopted a simplified version of the 


Besangon galaxy model (Robin et al. 2003 20041. A 


halo population at 100 times the measured density with 
{U, V, W) = (0.0, —220.0,0.0) kms“^ and velocity dispersion 
{au,av,iyw) = (131.0,106.0, 85.0) kms“^ was generated us¬ 
ing the same recipe described in Section The halo pop¬ 
ulation is then analyzed using the original set of kinematic 
properties. The same WDLF from Section [4.2| is overplotted 
for comparison in Fig. |10| The two WDLFs agree at all mag¬ 
nitudes and even the trends in underdensities at faint mag¬ 
nitudes are very similar. The differences in the {W/Vmax) 
distribution at faint magnitudes can be explained by the 
same argument as in Section |4.4[ Due to the similar kine¬ 
matic properties between the halo kinematic models applied 
in this work and the Besangon model, the size of this effect 
is much smaller. From this simple experiment we conclude 
that small differences between the assumed and underlying 
kinematics results in little differences in both the derived 
WDLF and the {V/Vmax) distribution. 


selected on the basis of both magnitude and proper mo¬ 
tion. Our simulations have shown that the assumption of 
setting the tangential velocity of an object as its intrinsic 
property would be invalid when a sample is drawn from a 
survey that has restricted proper motion limits and for pop¬ 
ulations where the tangential velocity distribution function 
varies over the sky. In order to select high quality proper 
motion objects when analysing real survey data the lower 
proper motion limit and tangential velocity limit are usu¬ 
ally set at such a level that the velocity space is very re¬ 
stricted. Therefore, it is necessary to correlate the tangential 
velocity with the proper motion depending on the proper¬ 
ties of each individual object. Two biases were identified 
from dissecting the generalized Ymax integral. The first one 
is an inconsistent density measurement of any non-uniform 
sample in which the density depends on the line of sight 
direction. The generalized method only considers the pop¬ 
ulation spatial density and survey limits, so different densi¬ 
ties would be measured at different line of sight directions 
as the distribution of proper motion and tangential velocity 
vary with line of sight. The other arises from the assump¬ 
tion that any given object would carry a constant tangential 
velocity. The choice of the upper proper motion limit has 
to be carefully determined from the properties of the sur¬ 
vey. An arbitrarily large upper limit has negligible effects 
in the analysis as the survey volume is sensitive only to the 
lower proper motion limit. The discovery fraction is, how¬ 
ever, sensitive to both upper and lower proper motion limits. 
From our simulations, it is found that surveys with upper 
proper motion limits smaller than ~0.5” yr“^ require more 
attention. This value seems small to the modern and future 
surveys (e.g. Pan-STARRS, Gaia, LSST etc.) but when a 
pairing is done with the early epoch Schmidt survey plate 
data for the bright objects for example the upper proper mo¬ 
tion limits would shrink rapidly when the maximum epoch 
differences can be up to more than half a century (the pair¬ 
ing radius cannot be increased indefinitely when associating 
catalogues with very different depth and resolution without 
introducing many spurious source matches). 

Population membership association has to be done care¬ 
fully because with a wrong set of kinematics the derived set 
of discovery fractions would become completely meaningless. 
However, in the case that a population completely domi¬ 
nates the density budget, for example the thin disc objects 
in the solar neighbourhood, the wrong membership asso¬ 
ciations have negligible effects to the luminosity functions. 
Finally, we would like to remind the readers to consider the 
limitations of (V/Ymax) carefully: (i) a fiat V/Vmax distri¬ 
bution with (V/Ymax) ~ 0.5 only indicates if the analysis 
is not biased towards distant or nearer objects, it gives no 
indication of the completeness of the sample. For a given 
flat distribution, one can conclude, at best, that the sur¬ 
vey is equally sensitive to objects at any distance. How¬ 
ever, (ii) given a complete sample, a flat distribution and 
(Y/Ymax) ~ 0.5 are expected. 


5 CONCLUSION 

We have presented an improved treatment of the proper 
motion and tangential velocity limits to arrive at an unbi¬ 
ased luminosity function estimator for stellar populations 


Application to other density estimators 

The treatment of detection limits due to proper motion se¬ 
lection have never been included in other density estima¬ 
tors. The modified Schmidt estimator is the only one that 
has been extended to incorporate such constraints. How- 
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ever, as detailed earlier, the treatment breaks down if we 
have restricted proper motion selection criteria. |Geijo et al.| 
120061 compared the properties of three different estima¬ 
tors: generalized Schmidt, Stepwise Maximum Likelihood 
and Choloniewski method. It is clear that the treatment of 
the proper motion limits were only applied to the Schmidt 
estimator by considering the tangential velocity as an in¬ 
trinsic property of an object, while the other two consid¬ 
ered the case of a magnitude limited sample. Their analysis 
may be affected by the modifications to the density estima¬ 
tor made in this work because the sample adopted in their 
work has assumed a fixed tangential velocity along line of 
sight. However, their work focused on the thin disc which 
is much less susceptible to the bias we identified, and by 
adopting a large upper proper motion limit in their data se¬ 
lection (/r £ [0.16, 2.00]), the effect should be tiny. Nonethe¬ 
less we would like to remind readers to pay attention to the 
faint end of their analysis which is most affected if the bias 
is noticeable at all. It should be possible to generalize the 
discovery fraction (as a completeness correction arising from 
tangential velocity and/or proper motion limit) to other den¬ 
sity estimators as the proper motion is now decoupled from 
the detection limit. 


Small Scale Variations 

Modern digital surveys are probing larger volumes 
with multi-epoch measurements over large areas of sky 
(e.g. SDSS, VST, Pan-STARRS, LSST and Gaia). Astro- 
metric precision is, however, a strong function of the direc¬ 
tion of the lines of sight due to varying observing qualities 
and imperfect camera fill factor arising from GCD chip gaps. 
In order to fully utilize the sample volume probed in any 
single or combined surveys, it is necessary to generalize the 
density estimators to cope with small scale fluctuations of 
the survey properties. Otherwise, the volume of the survey 
must be limited to the shallowest part of the sky. One pos¬ 
sible solution is to pixelise the sky into equal area pixels 
using HEALPi?!^ The photometric and astrometric preci¬ 
sion of each pixel are then limited to the worst observing 
condition within the pixel. 


Galactic model 

The current Schwarzschild distribution assumes a Gartesian 
system centred at the Sun. This approximation is good for 
surveys that probe only small distances (a few hundred par¬ 
secs). However, in the future surveys, most notably the Gaia 
and LSST, where the lower proper motion limits at the de¬ 
tection limits would be as low as 0.2masyr“^ [Goaia ~ 
20mag), and lmasyr“^ {tlsst ~ 24mag) respectively, an 
object with a tangential velocity of 20kms“^ could be de¬ 
tected with 5(7 confidence at ~17kpc and ~3.5kpc respec¬ 
tively. A WD with Tefi ~ 3000 K and Mi,oi ~ 16 mag, which 
translate to Goaia ~15.8mag and gnssx ~17.3mag at 10 pc, 
would be detectable at ~l-2 kpc in both cases. The overden¬ 
sity at the spiral arms is not accounted for in any work to 
date because the current surveys have not reached such dis¬ 
tances where the nearest spiral arm lies. However, when the 

^ http://healpix.jpl.nasa.gov/ 


faintest objects can be detected at a distance of over a thou¬ 
sand parsecs, the brighter objects would lie easily inside the 
spiral arm region which is ~lkpc from the Sun. In such a 
situation it is necessary to employ a more detailed density 
profile to account for the varying density in not just the ver¬ 
tical direction but also in the planar directions. With the 
decoupling of the proper motion limit from the photomet¬ 
ric limits, the incorporation of a more sophisticated model 
would simply be a translation of the Galactic density and 
velocity profiles from the galactocentric frame to the geo¬ 
centric frame. The key difference is that a Schwarzschild 
distribution function gives a constant tangential velocity 
distribution function as a function of the line of sight (ie. 
X = x(q^i^))i but in the adoption of a complex model, the 
density and velocity distribution functions have to be found 
as functions of both distance and the direction of line of 
sight such that x = tD). 
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